%% Graph: Quality Distribution

    fig = figure('visible','off','DefaultAxesFontSize',12);
    hold on
    scatter(const.Qgrid, eqm.nq, 15,'filled','o', 'MarkerFaceColor',[0.6,0.6,0.6])
    scatter(const.Qgrid, eqm_cf1.nq, 15,'filled','d', 'MarkerFaceColor',[0.3,0.3,0.3])
    hold off
    legend('Initial Equilibrium', 'Counterfactual', 'Location','northeast')
    xlabel('Quality')
    ylabel('Density')
    ax = gca;
    ax.YGrid = 'on';
    ax_pos = ax.Position;
    ax2 = axes('Position', ax_pos, 'XAxisLocation','top', 'YAxisLocation','right','Color','none','YTick',[]);
    xlim(ax2,[0, 8])
    xticks([0:8])
    %fprintf(' %.2f',const.nq_label_lnw)
    xticklabels({'','-0.25','0.04','0.28','0.49','0.68','0.91','1.06',''}) 
    xlabel('Log wage per worker (mean zero)')   
    %exportgraphics(fig,'../../output/figures/fullmodel_quality_distribution.pdf')    
    exportgraphics(fig,'../../output/figures/Figure5.pdf')    

    
     
%% Graph: Decomposition of domestic Pi(q) change    

    cost_component = (eqm_cf1.c./eqm.c).^(param.alpha_m*(1-param.sigma)*const.gamma)...
        .*(eqm_cf1.P_s/1).^(param.alpha_s*(1-param.sigma)*const.gamma)...
        .*(param.w/1).^((1-param.alpha_m-param.alpha_s)*(1-param.sigma)*const.gamma)...
        .*(param.w/1).^(1-const.gamma);
    demand_component = (eqm_cf1.D0./eqm.D0).^const.gamma;
    total = (eqm_cf1.Pi0./eqm.Pi0);
    check = total - cost_component.*demand_component;

    fig = figure('visible','off','DefaultAxesFontSize',12);
    hold on
    plot(const.Qgrid, cost_component,'LineStyle','-.','LineWidth',2.5, 'Color',ones(1,3)*0.3)
    plot(const.Qgrid, demand_component,'LineStyle','--','LineWidth',2.5, 'Color',ones(1,3)*0.4)
    plot(const.Qgrid, total,'LineWidth',3.5, 'Color',ones(1,3)*0.6)
    plot(const.Qgrid, ones(param.Q,1),'LineWidth',0.5,'Color','black')
    hold off
    %legend('C(1,q)^{\gamma(1-\sigma)}', 'D(q)^{\gamma}', '\Pi(q)', 'Location','northwest')
    legend('Cost Component', 'Demand Component', 'Total Domestic Profit Shifter', 'Location','northwest')
    ylabel('Counterfactual / Initial Equilibrium')
    xlabel('Quality')   
    ax = gca;
    ylim(ax,[0.9, 1.15]);  
    %exportgraphics(fig,'../../output/figures/fullmodel_Pi0_decomposition.pdf')    
    exportgraphics(fig,'../../output/figures/Figure6_A.pdf')    

    
    